#include "wrapper_hypre.h"
#include <vector>
typedef HYPRE_Int idx_t;
typedef HYPRE_Real data_t;

const idx_t num_diag = 5;

struct BlkInfo {
    int pid = -1;
    idx_t coords[2] = {-1, -1};
    idx_t begs[2] = {-1, -1};
    idx_t ends[2] = {-1, -1};
};

// 依照hypre的习惯，[0]在最内维，[1]在最外维
// const idx_t glb_dims [2] = { 3360, 3072 };
// const idx_t blk_dims [2] = {  160,  128 };
// const idx_t glb_dims [2] = { 11520, 14580 };
// const idx_t blk_dims [2] = {   960,   540 };

int main(int argc, char* argv[])
{
    int my_pid, num_procs;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_pid);
    MPI_Comm_size(MPI_COMM_WORLD, &num_procs);

    idx_t arg_cnt = 1;
    idx_t glb_dims[2], blk_dims[2];
    glb_dims[0] = atoi(argv[arg_cnt+1]);
    glb_dims[1] = atoi(argv[arg_cnt  ]);
    blk_dims[0] = atoi(argv[arg_cnt+3]);
    blk_dims[1] = atoi(argv[arg_cnt+2]); arg_cnt += 4;
    const std::string pathname(argv[arg_cnt++]);
    const std::string solver_name = std::string(argv[arg_cnt++]),
                        prec_name = std::string(argv[arg_cnt++]);

    assert(glb_dims[0] % blk_dims[0] == 0 && glb_dims[1] % blk_dims[1] == 0);
    std::vector<BlkInfo> blk_infos;

    const std::string distinfo_name = "distinfo." + std::to_string(glb_dims[0]) + "x" + std::to_string(glb_dims[1])
                                        + ".tile" + std::to_string(blk_dims[0]) + "x" + std::to_string(blk_dims[1]);
    FILE * fp = fopen((pathname + "/" + distinfo_name).c_str(), "r");
    fscanf(fp, "i1 i2 j1 j2 i j Index_withLand Index_rmvland");
    for (idx_t b = 0; b < (glb_dims[0] / blk_dims[0]) * (glb_dims[1] / blk_dims[1]); b++) {
        idx_t i1, i2, j1, j2, pi, pj, index_withLand, index_rmvLand;
        fscanf(fp, "%d,%d,%d,%d,%d,%d,%d,%d", & i1, & i2, & j1, & j2, & pi, & pj, & index_withLand, & index_rmvLand);
        if (index_rmvLand == -9999) continue;

        BlkInfo blk;
        blk.pid = blk_infos.size(); assert(blk.pid == index_rmvLand - 1);
        blk.coords[0] = pi - 1;
        blk.coords[1] = pj - 1;
        blk.begs[0] = i1 - 1;
        blk.ends[0] = i2; assert(blk.ends[0] <= glb_dims[0]);

        blk.begs[1] = j1 - 1;
        blk.ends[1] = j2; assert(blk.ends[1] <= glb_dims[1]);

        assert(blk.ends[0] - blk.begs[0] == blk_dims[0]);
        assert(blk.ends[1] - blk.begs[1] == blk_dims[1]);
        blk_infos.emplace_back(blk);
    }
    fclose(fp);
    
    const BlkInfo & my_blk = blk_infos[my_pid];
    const idx_t my_nelems = (my_blk.ends[0] - my_blk.begs[0]) * (my_blk.ends[1] - my_blk.begs[1]);
    // 先读入本进程所负责的区域的数据
    data_t * A_buf = new data_t [my_nelems * num_diag];
    data_t * rhs_buf = new data_t [my_nelems];// 读入数据时从外到内为 [3360][3072]
    data_t * x0_buf  = new data_t [my_nelems];// 读入数据时从外到内为 [3360][3072]

    // for (int p = 0; p < num_procs; p++) {
    //     if (p == my_pid) {
    //         printf("Proc %d : [ %d , %d ) x [ %d , %d )\n", my_pid, my_blk.begs[0], my_blk.ends[0], my_blk.begs[1], my_blk.ends[1]);
    //         fflush(stdout);
    //     }
    //     MPI_Barrier(MPI_COMM_WORLD);
    // }

    {// 读入向量
        assert(sizeof(double) == sizeof(HYPRE_Real));
        MPI_File fh = MPI_FILE_NULL;// 文件句柄
        MPI_Datatype etype = MPI_DOUBLE;
        MPI_Datatype read_type = MPI_DATATYPE_NULL;// 写出类型
        HYPRE_Int size[2], subsize[2], start[2];
        // 注意这里要反过来，因为存储的数据从外到内是[3072][3360]的
        size   [0] = glb_dims[1]   ; size   [1] = glb_dims[0];
        subsize[0] = blk_dims[1]   ; subsize[1] = blk_dims[0];
        start  [0] = my_blk.begs[1]; start  [1] = my_blk.begs[0];

        MPI_Type_create_subarray(2, size, subsize, start, MPI_ORDER_C, etype, &read_type);
        MPI_Type_commit(&read_type);
        MPI_Offset displacement = 0;
        displacement *= sizeof(HYPRE_Real);// 位移要以字节为单位！
        MPI_Status status;

        int ret;
        ret = MPI_File_open(MPI_COMM_WORLD, (pathname + "/array_b").c_str(), MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);
        if (ret != MPI_SUCCESS) {
            printf("Could not open b file: ret %d\n", ret);
        }
        MPI_File_set_view(fh, displacement, etype, read_type, "native", MPI_INFO_NULL);
        MPI_File_read_all(fh, rhs_buf, my_nelems, etype, &status);
        MPI_File_close(&fh);

        ret = MPI_File_open(MPI_COMM_WORLD, (pathname + "/array_x0").c_str(), MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);
        if (ret != MPI_SUCCESS) {
            printf("Could not open x file: ret %d\n", ret);
        }
        MPI_File_set_view(fh, displacement, etype, read_type, "native", MPI_INFO_NULL);
        MPI_File_read_all(fh, x0_buf, my_nelems, etype, &status);
        MPI_File_close(&fh);

        MPI_Type_free(&read_type);
    }
    {// 读入矩阵
        assert(sizeof(double) == sizeof(HYPRE_Real));
        MPI_File fh = MPI_FILE_NULL;// 文件句柄
        MPI_Datatype etype = MPI_DOUBLE;
        MPI_Datatype read_type = MPI_DATATYPE_NULL;// 写出类型
        HYPRE_Int size[3], subsize[3], start[3];
        // 注意这里要反过来，因为存储的数据从外到内是[3072][3360]的
        size   [0] = glb_dims[1]   ; size   [1] = glb_dims[0];    size   [2] = num_diag;
        subsize[0] = blk_dims[1]   ; subsize[1] = blk_dims[0];    subsize[2] = num_diag;
        start  [0] = my_blk.begs[1]; start  [1] = my_blk.begs[0]; start  [2] = 0;

        MPI_Type_create_subarray(3, size, subsize, start, MPI_ORDER_C, etype, &read_type);
        MPI_Type_commit(&read_type);
        MPI_Offset displacement = 0;
        displacement *= sizeof(HYPRE_Real);// 位移要以字节为单位！
        MPI_Status status;

        int ret;
        ret = MPI_File_open(MPI_COMM_WORLD, (pathname + "/mat.AOS.bin").c_str(), MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);
        if (ret != MPI_SUCCESS) {
            printf("Could not open file: ret %d\n", ret);
        }
        MPI_File_set_view(fh, displacement, etype, read_type, "native", MPI_INFO_NULL);
        MPI_File_read_all(fh, A_buf, my_nelems * num_diag, etype, &status);
        MPI_File_close(&fh);
        MPI_Type_free(&read_type);
    }

    // 手动检验
    {
        double loc_b_dot = 0.0, loc_x_dot = 0.0, loc_A_dot = 0.0;
        for (idx_t i = 0; i < my_nelems; i++) {
            loc_b_dot += rhs_buf[i] * rhs_buf[i];
            loc_x_dot +=  x0_buf[i] *  x0_buf[i];
            if (A_buf[i * num_diag + 2] == -1.0) continue;
            // for (idx_t d = 0; d < num_diag; d++)
                loc_A_dot += A_buf[i * num_diag + 2] * A_buf[i * num_diag + 2];
        }
        // printf("Proc %d : %.15e %.15e\n", my_pid, loc_b_dot, loc_x_dot);
        double ck_b_dot = 0.0, ck_x_dot = 0.0, ck_A_dot = 0.0;
        MPI_Allreduce(& loc_b_dot, & ck_b_dot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        MPI_Allreduce(& loc_x_dot, & ck_x_dot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        MPI_Allreduce(& loc_A_dot, & ck_A_dot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        if (my_pid == 0) {
            printf("My Calc:\n");
            printf("(  b,   b) = %.20e\n",  (double)ck_b_dot);
            printf("(  x,   x) = %.20e\n",  (double)ck_x_dot);
            printf("(  A,   A) = %.20e\n",  (double)ck_A_dot);
        }
    }

    idx_t obj_type = HYPRE_SSTRUCT;
    if (prec_name == "AMG") obj_type = HYPRE_PARCSR;
    else if (prec_name == "SMG" || prec_name == "PFMG") obj_type = HYPRE_STRUCT;
    if (my_pid == 0) {
        printf("object type : %d\n", obj_type);
    }

    HYPRE_Init();

    HYPRE_SStructGrid ssgrid;
    HYPRE_SStructStencil stencil;
    HYPRE_SStructGraph ssgraph;
    HYPRE_SStructMatrix   A;
    HYPRE_SStructVector   b, x, y, x0;
    HYPRE_ParCSRMatrix par_A;
    HYPRE_ParVector par_b, par_x, par_y;
    HYPRE_StructMatrix s_A;
    HYPRE_StructVector s_b, s_x, s_y;
    TEST_RECORD records[TEST_CNT];
    TEST_CONFIG config;
    config.case_name = "ocean";
    config.rtol = 1.0e-7;

    {
        idx_t nparts = -1;
        idx_t my_part = -1;
        idx_t my_ilower[2] = {-1, -1}, my_iupper[2] = {-1, -1};
        if (prec_name == "PFMG" || prec_name == "SMG") {
            nparts = 1;
            my_part = 0;
            if (my_pid == 0) printf("  Total #Blks : %d\n", nparts);

            HYPRE_SStructGridCreate(MPI_COMM_WORLD, 2, nparts, &ssgrid);
            my_ilower[0] = my_blk.begs[0]; my_iupper[0] = my_blk.ends[0] - 1;
            my_ilower[1] = my_blk.begs[1]; my_iupper[1] = my_blk.ends[1] - 1;// 注意是左闭右闭区间
            HYPRE_SStructGridSetExtents(ssgrid, 0, my_ilower, my_iupper);

            HYPRE_SStructVariable vtypes[1] = { HYPRE_SSTRUCT_VARIABLE_CELL };
            for (idx_t part = 0; part < nparts; part++) {// 每个part都要执行定义变量的操作
                HYPRE_SStructGridSetVariables(ssgrid, part, sizeof(vtypes)/sizeof(HYPRE_SStructVariable), vtypes);
            }
            // 单块无非stencil关系
        }
        else {
            nparts = blk_infos.size();
            my_part = my_blk.pid;
            assert(nparts == num_procs);
            if (my_pid == 0) printf("  Total #Blks : %d\n", nparts);

            HYPRE_SStructGridCreate(MPI_COMM_WORLD, 2, nparts, &ssgrid);
            my_ilower[0] = 0; my_iupper[0] = my_blk.ends[0] - my_blk.begs[0] - 1;
            my_ilower[1] = 0; my_iupper[1] = my_blk.ends[1] - my_blk.begs[1] - 1;// 注意是左闭右闭区间
            HYPRE_SStructGridSetExtents(ssgrid, my_part, my_ilower, my_iupper);

            HYPRE_SStructVariable vtypes[1] = { HYPRE_SSTRUCT_VARIABLE_CELL };
            for (idx_t part = 0; part < nparts; part++) {// 每个part都要执行定义变量的操作
                HYPRE_SStructGridSetVariables(ssgrid, part, sizeof(vtypes)/sizeof(HYPRE_SStructVariable), vtypes);
            }
            // 添加上非stencil的连接关系
            {
            idx_t index_maps[2] = {0, 1};
            idx_t index_dirs[2] = {1, 1};

            // 第0方向的邻居 {0, -1}
            for (size_t k = 0; k < blk_infos.size(); k++) {
                const BlkInfo & ngb_blk = blk_infos[k];
                if (ngb_blk.begs[0] == my_blk.begs[0] &&// 相同
                    ngb_blk.ends[0] == my_blk.ends[0] &&// 相同
                    ngb_blk.ends[1] == my_blk.begs[1] ) {// 连续
                    // assert(ngb_blk.coords[0] == my_blk.coords[0]);
                    // assert(ngb_blk.coords[1] == my_blk.coords[1] - 1);
                    idx_t halo_ilower[2] = { my_ilower[0], my_ilower[1] - 1 };
                    idx_t halo_iupper[2] = { my_iupper[0], my_ilower[1] - 1 };
                    idx_t ngb_ilower [2] = { my_ilower[0], my_iupper[1] };
                    idx_t ngb_iupper [2] = { my_iupper[0], my_iupper[1] };
#ifdef _zyDEBUG
                    // printf("Proc %d: part %d (%d,%d)~(%d,%d) with ngb_part %d (%d,%d)~(%d,%d)\n",
                    //     my_pid, my_blk.pid, 
                    //     halo_ilower[0], halo_ilower[1],
                    //     halo_iupper[0], halo_iupper[1],
                    //     k,
                    //     ngb_ilower[0], ngb_ilower[1],
                    //     ngb_iupper[0], ngb_iupper[1]
                    // );
#endif
                    HYPRE_SStructGridSetNeighborPart(ssgrid, my_blk.pid, halo_ilower, halo_iupper,
                            k, ngb_ilower, ngb_iupper, index_maps, index_dirs);   
                    break;// 找到该方向的邻居后不必再找   
                }
            }
            // 第1方向的邻居 {-1, 0}
            for (size_t k = 0; k < blk_infos.size(); k++) {
                const BlkInfo & ngb_blk = blk_infos[k];
                if (ngb_blk.begs[1] == my_blk.begs[1] &&// 相同
                    ngb_blk.ends[1] == my_blk.ends[1] &&// 相同
                    ngb_blk.ends[0] == my_blk.begs[0] ) {// 连续
                    // assert(ngb_blk.coords[1] == my_blk.coords[1]);
                    // assert(ngb_blk.coords[0] == my_blk.coords[0] - 1);
                    idx_t halo_ilower[2] = { my_ilower[0] - 1, my_ilower[1] };
                    idx_t halo_iupper[2] = { my_ilower[0] - 1, my_iupper[1] };
                    idx_t ngb_ilower [2] = { my_iupper[0]    , my_ilower[1] };
                    idx_t ngb_iupper [2] = { my_iupper[0]    , my_iupper[1] };
#ifdef _zyDEBUG
                    // printf("Proc %d: part %d (%d,%d)~(%d,%d) with ngb_part %d (%d,%d)~(%d,%d)\n",
                    //     my_pid, my_blk.pid, 
                    //     halo_ilower[0], halo_ilower[1],
                    //     halo_iupper[0], halo_iupper[1],
                    //     k,
                    //     ngb_ilower[0], ngb_ilower[1],
                    //     ngb_iupper[0], ngb_iupper[1]
                    // );
#endif
                    HYPRE_SStructGridSetNeighborPart(ssgrid, my_blk.pid, halo_ilower, halo_iupper,
                            k, ngb_ilower, ngb_iupper, index_maps, index_dirs);   
                    break;// 找到该方向的邻居后不必再找   
                }
            }
            // 第3方向的邻居 { 1, 0}
            for (size_t k = 0; k < blk_infos.size(); k++) {
                const BlkInfo & ngb_blk = blk_infos[k];
                if (ngb_blk.begs[1] == my_blk.begs[1] &&// 相同
                    ngb_blk.ends[1] == my_blk.ends[1] &&// 相同
                    ngb_blk.begs[0] == my_blk.ends[0] ) {// 连续
                    // assert(ngb_blk.coords[1] == my_blk.coords[1]);
                    // assert(ngb_blk.coords[0] == my_blk.coords[0] + 1);
                    idx_t halo_ilower[2] = { my_iupper[0] + 1, my_ilower[1] };
                    idx_t halo_iupper[2] = { my_iupper[0] + 1, my_iupper[1] };
                    idx_t ngb_ilower [2] = { my_ilower[0]    , my_ilower[1] };
                    idx_t ngb_iupper [2] = { my_ilower[0]    , my_iupper[1] };
#ifdef _zyDEBUG
                    // printf("Proc %d: part %d (%d,%d)~(%d,%d) with ngb_part %d (%d,%d)~(%d,%d)\n",
                    //     my_pid, my_blk.pid, 
                    //     halo_ilower[0], halo_ilower[1],
                    //     halo_iupper[0], halo_iupper[1],
                    //     k,
                    //     ngb_ilower[0], ngb_ilower[1],
                    //     ngb_iupper[0], ngb_iupper[1]
                    // );
#endif
                    HYPRE_SStructGridSetNeighborPart(ssgrid, my_blk.pid, halo_ilower, halo_iupper,
                            k, ngb_ilower, ngb_iupper, index_maps, index_dirs);   
                    break;// 找到该方向的邻居后不必再找   
                }
            }
            // 第4方向的邻居 { 0, 1}
            for (size_t k = 0; k < blk_infos.size(); k++) {
                const BlkInfo & ngb_blk = blk_infos[k];
                if (ngb_blk.begs[0] == my_blk.begs[0] &&// 相同
                    ngb_blk.ends[0] == my_blk.ends[0] &&// 相同
                    ngb_blk.begs[1] == my_blk.ends[1] ) {// 连续
                    // assert(ngb_blk.coords[0] == my_blk.coords[0]);
                    // assert(ngb_blk.coords[1] == my_blk.coords[1] + 1);
                    idx_t halo_ilower[2] = { my_ilower[0], my_iupper[1] + 1 };
                    idx_t halo_iupper[2] = { my_iupper[0], my_iupper[1] + 1 };
                    idx_t ngb_ilower [2] = { my_ilower[0], my_ilower[1] };
                    idx_t ngb_iupper [2] = { my_iupper[0], my_ilower[1] };
#ifdef _zyDEBUG
                    // printf("Proc %d: part %d (%d,%d)~(%d,%d) with ngb_part %d (%d,%d)~(%d,%d)\n",
                    //     my_pid, my_blk.pid, 
                    //     halo_ilower[0], halo_ilower[1],
                    //     halo_iupper[0], halo_iupper[1],
                    //     k,
                    //     ngb_ilower[0], ngb_ilower[1],
                    //     ngb_iupper[0], ngb_iupper[1]
                    // );
#endif
                    HYPRE_SStructGridSetNeighborPart(ssgrid, my_blk.pid, halo_ilower, halo_iupper,
                            k, ngb_ilower, ngb_iupper, index_maps, index_dirs);   
                    break;// 找到该方向的邻居后不必再找   
                }
            }
        }
        }
        HYPRE_SStructGridAssemble(ssgrid);// a collective call finalizing the grid assembly.

        idx_t offsets[num_diag][2] = {
            {0, -1}, {-1, 0}, {0,0}, {1, 0}, {0, 1}
            // {-1, 0}, {0, -1}, {0,0}, {0, 1}, {1, 0}
        };
        HYPRE_SStructStencilCreate(2, num_diag, &stencil);
        for (idx_t e = 0; e < num_diag; e++) {
            HYPRE_SStructStencilSetEntry(stencil, e, offsets[e], 0);// 只有0号变量一种
        }

        HYPRE_SStructGraphCreate(MPI_COMM_WORLD, ssgrid, &ssgraph);
        HYPRE_SStructGraphSetObjectType(ssgraph, obj_type);
        for (idx_t part = 0; part < nparts; part++) {// 每个part都要执行定义模板的操作
            HYPRE_SStructGraphSetStencil(ssgraph, part, 0, stencil);
        }
        HYPRE_SStructGraphAssemble(ssgraph);

        // 建立向量
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &b);// Create an empty vector object
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &x);
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &y);
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &x0);
        HYPRE_SStructVectorSetObjectType(b, obj_type);// Set the object type for the vectors to be the same as was already set for the matrix
        HYPRE_SStructVectorSetObjectType(x, obj_type);
        HYPRE_SStructVectorSetObjectType(y, obj_type);
        HYPRE_SStructVectorSetObjectType(x0, obj_type);
        HYPRE_SStructVectorInitialize(b);// Indicate that the vector coefficients are ready to be set
        HYPRE_SStructVectorInitialize(x);
        HYPRE_SStructVectorInitialize(y);
        HYPRE_SStructVectorInitialize(x0);
        {// 填充向量数据
            // HYPRE_SStructVectorSetBoxValues(x, my_part, my_ilower, my_iupper, 0, x0_buf);
            HYPRE_SStructVectorSetConstantValues(x, 0.0);
            HYPRE_SStructVectorSetBoxValues(b, my_part, my_ilower, my_iupper, 0, rhs_buf);

            HYPRE_SStructVectorSetConstantValues(y, 0.0);
        }
        HYPRE_SStructVectorAssemble(b);
        HYPRE_SStructVectorAssemble(x);
        HYPRE_SStructVectorAssemble(y);// 对y也要创建
        HYPRE_SStructVectorAssemble(x0);
        HYPRE_SStructVectorCopy(x, x0);

        HYPRE_SStructMatrixCreate(MPI_COMM_WORLD, ssgraph, &A);// Create an empty matrix object
        HYPRE_SStructMatrixSetObjectType(A, obj_type);
        HYPRE_SStructMatrixInitialize(A);// Indicate that the matrix coefficients are ready to be set 
        {// 填充矩阵数据
            idx_t stencil_indices[num_diag];
            for (idx_t j = 0; j < num_diag; j++)
                stencil_indices[j] = j;
            HYPRE_SStructMatrixSetBoxValues(A, my_part, my_ilower, my_iupper, 0, num_diag, stencil_indices, A_buf);
        }
        HYPRE_SStructMatrixAssemble(A);// a collective call finalizing the matrix assembly.

        for (idx_t test = 0; test < TEST_CNT; test++) {
            HYPRE_SStructVectorCopy(x0, x);

            data_t b_dot, x_dot, Ab_dot, Ax_dot;
            idx_t ret;
            // 做spmv检验一下数据传对没有
            if (obj_type == HYPRE_PARCSR) {
                HYPRE_SStructMatrixGetObject(A, (void **) &par_A);
                HYPRE_SStructVectorGetObject(b, (void **) &par_b);
                HYPRE_SStructVectorGetObject(x, (void **) &par_x);
                HYPRE_SStructVectorGetObject(y, (void **) &par_y);
                ret = HYPRE_ParVectorInnerProd(par_b, par_b, &b_dot);
                ret = HYPRE_ParVectorInnerProd(par_x, par_x, &x_dot);
                HYPRE_ParCSRMatrixMatvec(1.0, par_A, par_b, 0.0, par_y);
                ret = HYPRE_ParVectorInnerProd(par_y, par_y, &Ab_dot);
                HYPRE_ParCSRMatrixMatvec(1.0, par_A, par_x, 0.0, par_y);
                ret = HYPRE_ParVectorInnerProd(par_y, par_y, &Ax_dot);
            } else if (obj_type == HYPRE_STRUCT) {
                HYPRE_SStructMatrixGetObject(A, (void **) &s_A);
                HYPRE_SStructVectorGetObject(b, (void **) &s_b);
                HYPRE_SStructVectorGetObject(x, (void **) &s_x);
                HYPRE_SStructVectorGetObject(y, (void **) &s_y);

                HYPRE_StructVector s_x0;
                HYPRE_SStructVectorGetObject(x0, (void **) &s_x0);
                HYPRE_StructVectorCopy(s_x0, s_x);

                ret = HYPRE_StructVectorInnerProd(s_b, s_b, &b_dot);
                ret = HYPRE_StructVectorInnerProd(s_x, s_x, &x_dot);
                HYPRE_StructMatrixMatvec(1.0, s_A, s_b, 0.0, s_y);
                ret = HYPRE_StructVectorInnerProd(s_y, s_y, &Ab_dot);
                HYPRE_StructMatrixMatvec(1.0, s_A, s_x, 0.0, s_y);
                ret = HYPRE_StructVectorInnerProd(s_y, s_y, &Ax_dot);
            } else { assert(obj_type == HYPRE_SSTRUCT);
                ret = HYPRE_SStructInnerProd(b, b, &b_dot);
                ret = HYPRE_SStructInnerProd(x, x, &x_dot);
                HYPRE_SStructMatrixMatvec(1.0, A, b, 0.0, y);
                ret = HYPRE_SStructInnerProd(y, y, &Ab_dot);
                HYPRE_SStructMatrixMatvec(1.0, A, x, 0.0, y);
                ret = HYPRE_SStructInnerProd(y, y, &Ax_dot);
            }

            if (my_pid == 0) {
                printf("(  b,   b) = %.20e\n",  (double)b_dot);
                printf("(  x,   x) = %.20e\n",  (double)x_dot);
                printf("(A*b, A*b) = %.20e\n", (double)Ab_dot);
                printf("(A*x, A*x) = %.20e\n", (double)Ax_dot);
            }

            buildup_solver(solver_name, prec_name, config);

            idx_t num_iterations = 0;
            data_t t_setup = 0.0, t_solve = 0.0, final_res_norm = 0.0;

            if (obj_type == HYPRE_SSTRUCT)
                setup_and_solve(solver_name, prec_name, A, b, x, t_setup, t_solve, final_res_norm, num_iterations);
            else if (obj_type == HYPRE_STRUCT) {
                setup_and_solve(solver_name, prec_name, s_A, s_b, s_x, t_setup, t_solve, final_res_norm, num_iterations);
            }
            else if (obj_type == HYPRE_PARCSR) {
                setup_and_solve(solver_name, prec_name, par_A, par_b, par_x, t_setup, t_solve, final_res_norm, num_iterations);
                HYPRE_SStructVectorGather(b);
                HYPRE_SStructVectorGather(x);
            }

            // 计算真实残差
            data_t r_nrm2 = 0.0, b_nrm2 = 0.0;
            check_residual(prec_name, A, x, b, y, r_nrm2, b_nrm2);
            if (my_pid == 0) {
                printf("\033[1;35mtrue ||r|| = %20.16e ||r||/||b||= %20.16e\033[0m\n", r_nrm2, r_nrm2/b_nrm2);
                printf("Iterations = %d\n", num_iterations);
                printf("Time cost %.5f %.5f %.5f %d\n", (double)t_setup, (double)t_solve, (double)(t_setup + t_solve), num_iterations);
                printf("Final Relative Residual Norm = %e\n", (double)final_res_norm);
                printf("\n");
            }

            records[test].iter = num_iterations;
            records[test].setup = t_setup;
            records[test].solve = t_solve;

            destroy_solver(solver_name, prec_name);
        }

        if (my_pid == 0) {
            std::sort(records, records + TEST_CNT);
            const TEST_RECORD & best = records[0], & worst = records[TEST_CNT - 1];
            TEST_RECORD avg;
            for (idx_t i = 0; i < TEST_CNT; i++) {
                avg.setup += records[i].setup;
                avg.solve += records[i].solve;
                avg.prec  += records[i].prec;
                avg.iter  += records[i].iter;
            }
            avg.setup /= TEST_CNT;
            avg.solve /= TEST_CNT;
            avg.prec  /= TEST_CNT;
            avg.iter  /= TEST_CNT;
            data_t min_tot = best .setup + best .solve;
            data_t max_tot = worst.setup + worst.solve;
            data_t avg_tot = avg  .setup + avg  .solve;
            printf("\n  Summary of %d tests:\n", TEST_CNT);
            printf("     Setup time  Solve time  #Iter  Total time\n");
            printf("Min    %.6f    %.6f  %5d  %.6f\n", best .setup, best .solve, best .iter, min_tot);
            printf("Avg    %.6f    %.6f  %5d  %.6f\n", avg  .setup, avg  .solve, avg  .iter, avg_tot);
            printf("Max    %.6f    %.6f  %5d  %.6f\n", worst.setup, worst.solve, worst.iter, max_tot);    
        }

        HYPRE_SStructMatrixDestroy(A);
        HYPRE_SStructVectorDestroy(b); HYPRE_SStructVectorDestroy(x); HYPRE_SStructVectorDestroy(y); HYPRE_SStructVectorDestroy(x0);
        HYPRE_SStructGraphDestroy(ssgraph);
        HYPRE_SStructStencilDestroy(stencil);
        HYPRE_SStructGridDestroy(ssgrid);
    }
    HYPRE_Finalize();

    delete [] A_buf; delete [] rhs_buf; delete [] x0_buf;
    MPI_Finalize();
    return 0;
}